THE STATISTICAL MECHANICS OF THE SELF-GRAVITATING GAS: 
EQUATION OF STATE AND FRACTAL DIMENSION 
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We provide a complete picture of the self-gravitating non-relativistic gas at thermal equilibrium 
using Monte Carlo simulations (MC), analytic mean field methods (MF) and low density expansions. 
' The system is shown to possess an infinite volume limit, both in the canonical (CE) and in the mi- 

f*"^ , crocanonical ensemble (MCE) when N, V — » oo, keeping N/V 1/l3 fixed. We compute the equation 

■ of state (we do not assume it as is customary) , the entropy, the free energy, the chemical potential, 

£Nj ' the specific heats, the compressibilities, the speed of sound and analyze their properties, signs and 

singularities. The MF equation of state obeys a first order non-linear differential equation of Abel 
type. The MF gives an accurate picture in agreement with the MC simulations both in the CE and 
MCE. The inhomogeneous particle distribution in the ground state suggest a fractal distribution 
' with Haussdorf dimension D with D slowly decreasing with increasing density, 1 < D < 3. 

I> ; 

The study of the self-gravitating thermal gas has both fundamental and practical physical interest ||, |): (i) It 
possess remarkable thermodynamic properties due to the long range of the gravitational force [such properties never 
' occur in ordinary systems with short range forces], (ii) It plays a central role in astrophysics and cosmology; cold 
clouds in the interstellar medium and its remarkable observed scaling laws (from lCP 4 pc up to lOOpc), as well as 
galaxy distributions (up to 100 — 200Mpc) can be described with it. (hi) It is relevant for the study of stellar objects 
in the non-relativistic and relativistic cases. 

The ground state is inhomogeneous and the usual thermodynamic limit: number of particles N — > oo, volume 
V — > oo, with N/V fixed leads to collapse into a very dense phase. The gaseous phase can only exist when N, V — > oo 
Q\ . with N/V 1 / 3 fixed. This is a diluted limit where the particle density N/V goes to zero as V~ 2 I 3 . The appropriate 
dimensionless variable (for particles of mass m at temperature T on a box of volume V) is r) = G ™ T N (it can be a 
spherical or cubic box with volume V — L 3 or V — (4tt/3)R 3 , respectively). r\ is related to the Jeans' length of the 
O ,' gas dj through r\ = 3(L/dj) 2 . We call 'thermodynamic limit' the limit TV — > oo, R — > oo with a fixed ratio R/N. In 
D . this limit, physical magnitudes are expressed naturally as functions of r\. 
*^ For small r\, the gas behaves as a perfect gas. For growing 77, PV/NT (where P is the pressure) decreases (see 

, figs. 1-2) due to the attractive character of gravity. Finally, at some critical r\ C rit the gas exhibits a sharp clumping 
transition to a dense phase with negative pressure. The extension of the gaseous phase and the value r\ cr it depend 
on the thermodynamical ensemble (see figs. 1-2): the gas phase is larger in the microcanonical ensemble (MCE) and 
smaller in the canonical ensemble (CE). We investigate the self-gravitating gas with Monte Carlo (MC) and analytic 
mean field (MF) methods; in the dilute limit we expand in powers of ?/. Our results show that the CE and MCE 
yield the same results in their common range of the gaseous phase. The MF correctly describes the thermodynamic 
limit except near the critical points, the MF is valid for N\r\ — r\c.rit\ 3> 1- The vicinity of the critical point should be 
studied in a double scaling limit N — > 00, rj — * r) C rit- 

The particle distribution p(q) proves to be inhomogeneous (except for rj -C 1) and described by an universal 
function of the geometry, 77 and the ratio r = q/R. D slowly decreases from the value D = 3 for the ideal gas (77 = 0) 
till D = 0.98 in the extreme limit of the MC point taking the value 1.6 at r\ cr it [see Tabic 1]. The particle density in 
the bulk behaves as p(q) ~ r D ~ 3 . This indicates the presence of a fractal distribution with Haussdorf dimension D. 
In the MCE, the entropy can be written as 
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where E is the total energy, G is Newton's gravitational constant and 
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In the CE the partition function can be written as 
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We make now explicit the volume dependence by introducing the variables f%, l<l<Nasqi=Rf f i, ?i = 
(xi,yi,zi), where |r;| < 1 for a spherical box and < Xi,yi,zi < 1 for a cubic box. The momentum integrals can be 
straightforwardly computed both in eq. ([!]) and in eq. (Q) yielding 
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- — are dimensionless and we introduced the notation [X] 7 _ 



where u(fi, ...,r N ) = jj J2i<k 3 <n pp^ and £ = Gm^jv^ 
X" 9(X), being 9(X) the step function. All space integrals in eq.(||) and below are over an unit volume. The 
temperature is given in the MCE as a function of E and £ by (for N ^> 1) 
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In the CE, the free energy results 
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where F = — iVTlog ( 1 ^-) 3 ^ 2 is the free energy for an ideal gas. Defining f(rf) = 1 — |^ $' N (rj), we then find 
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P stands for the external pressure on the gas. Notice that the equation of state PV — NTf(rj) is here computed 
from the gravitational interaction, we do not assume an equation of state as it is customary. It follows from eqs.(|J) 
and (|) that ??£ = g(£) = 3 [f(rj) - 1/2] and that the virial theorem, (3PV = |iVT + E), holds both in the MCE and 
in the CE. 

We find that the following properties hold in the N — > oo , L — > oo limit provided 77 and £ are kept fixed: 

(a) f(j]) and g(£) reach finite limits. [In fact, fixed rj implies a finite limit for £ and viceversa]. 

(b) The energy, free energy, entropy and PV are proportional to N (besides a log N term in the entropy as for an 
ideal gas). 

(c) Both ensembles, the MCE and CE yield the same physical magnitudes. 

We have verified properties (a)-(c) in three ways. First, by direct calculation in the dilute regime (£ 3> 1 , T] <C 1) 
by expanding eqs.(||) and eqs.(||) in powers of l/£ and 77 for the MCE and CE, respectively. Second, by performing 
Monte Carlo simulations both in the MCE and the CE. Third, by mean field approximations to eqs.(|j) and eqs.(||). 

The specific heat at constant volume takes the form, 



CV = N{df) v = 3 
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This expression is valid in the thermodynamic limit both in the CE and MCE. In the CE, the specific heat is related 
to the fluctuations of the potential energy (AC/) 2 and it is positive defined, 



(cv)ce = ^ + (AC/) 2 , (AC/) 2 
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Dilute Limit Calculations. 

For a dilute gas, we expand the integrands in eqs.(|4|) and eqs.(^J) in powers of l/£ or 77, respectively. In the N — > 00 
limit the expressions simplify considerably. Moreover, divergent integrals in the zero cutoff limit are eliminated by 
1/N factors. That is, after the thermodynamic limit is taken the cutoff corrections are of the order 0(a 2 ) and clearly 
vanish in the a — > limit. We find after calculation, 
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For a sphere of unit volume we find b s ph = 1/5 , bf h = 153/35, whereas for a unit volume cube b c ube = 0.19462 .... 
Notice that the cube and sphere values differ only by about 3%. We see from eq.(||) that the perfect gas behaviour 
gets corrected by negative definite terms due to the attractive nature of the gravitational interaction. [The full 
correction to the ideal gas, $^(77), is positive definite according to eq.(@)]. 



Monte Carlo Calculations. 



We applied the standard Metropolis algorithm in a cube of size L with total energy E in the MCE and at temperature T 
in the CE. The number of particles N went up to 2000[^. We introduced a small short distance cutoff ~ 10 _3 i— 10~ 6 £ 
in the attractive Newton's potential. All results in the gaseous phase were insensitive to the cutoff value. 

Two different phases show up: a non-perfect gas for rj < rj c , and a condensed system with negative pressure for 
rj > r\ c . The transition between the two phases is very sharp, with a negative jump in the entropy from the gas to the 
condensed phase. This phase transition is associated with the Jeans instability. 

We plot in figs. 1-2 /(ry) = PV/[NT] as a function of rj. For small 7/, the MC results for PV/[NT] reproduce very 
well the analytical formula (||). PV/[NT] monotonically decreases with 77 as forecasted by the dilute expansion (||). 
For r\ = rjT — 1.51 (point T in fig. 2), close to rjc — 1-54 (point C in fig. 2) a phase transition suddenly happens in the 
CE and PV/[NT) becomes large and negative. The interparticle distance < r > monotonically decreases with rj too. 
When rj crosses rjT , < r > has a sharp decrease. 

The MCE and CE Monte Carlo results are very close (up to the statistical error) for < ij < rjc, that is for 
00 > £ > £c — —0.19. In the MCE the gas does not clump at rj = rjc (point C in fig. 2) and the specific heat becomes 
negative between the points C and MC. The gas does clump in the MCE at £ ~ —0.32 , 77 ~ 1.35 (point MC in fig. 2) 
increasing both its temperature and pressure discontinuously. As is clear, the domain between C and MC cannot be 
reached in the CE since cy > in the CE. 

Mean Field Calculations. 

We now recast the coordinate partition function e*™^ as a functional integral in the thermodynamic limit. 

e <M77)^l J J Dpdae -Ns[p(.)] + l a($ d 3 xp(3)-l) ^ 

s[p{.)\ = -\ f^^r P {x)p{y) + [d 3 xp(x)log[p(x)/e] . 
2 J \x - y\ J 

where we used the coordinates x in the unit volume. The first term is the potential energy, the second term is the 
functional integration measure for this case (see ||). Eq.(||) is dominated for large N by its stationary point solutions: 
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a is a Lagrange multiplier enforcing the constraint J d 3 x p(x) — 1. Applying the Laplacian and setting cb(x) = logp(x) 
yields, for the spherically symmetric case 
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with boundary conditions cb'(0) = and <f>'(l) = —J] 1 
defined as r, R = = r) (f) 1/3 = 1.61199 ... 77 

Using the scale covariance of eq.(nl|) M 



Here the variable 77^ appropriate for a spherical symmetry is 



) [|2j, 4>{r) can be expressed as cb(r) = log(A 2 /4 ir r/ R ) + x(Ar) where 
x"(A) + |x'(A)+e* (A) =0, x'(0)-0 



This equation is invariant under the transformation: 

A => A e Q 
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where a is a real number. Hence, we can set x(0) = without loosing generality. 



x(x) is independent of 77" and A is related to 77^ by A x'(A) 
Evaluating the functional integral in eq.(|^) by saddle point yields 
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where Dc(r] R ) stands for the determinant of small fluctuations around the spherically symmetric saddle point ( p"l| ) 
and C(rj R ) for the two-loop corrections. Dc(i] R ) can be expressed as an infinite product over the partial waves 
D c {r] R ) = n i>0 [Di(ri R )] 2l+1 . D {ii R ) and -Di^) can be computed in closed form g 

D 0(V R ) = A[A 2 (^)e x(A(r,ii)) - V R ] , D 1 ( V R ) = e*^ R )) . 
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Do(rj R ) [see fig. 1] is positive for < r) R < T] R = 2.51755 . . . and vanish linearly in J rj R — rj R at 77^ = rj R . 
We thus get for the free energy in the CE from eqs.(||), ([))), (10) and (|lj) 
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/MF(7? «) = ^ e *(^ B )) 

s{v R ) = 3[1 - f MF (v R )} -V R + \og[3f MF ( V R )/4n] 
It follows from cq.(|T^) and ( [l6| ) that /a/f('7' R ) obeys the first order equation 

^(3/-l)/'(77 i? ) + (3/-3 + 77 fi )/ = 0. 

which reduces to an Abel equation of first kindQ. 

By expanding the solution of eq.(|l^) in powers of A one checks that eq. (||) reproduces eq.(|){§. We plot in fig.l 
Jmf{ii R ) as a function of 77^ obtained by solving eqs.(]l^-|l5|) by the Runge-Kutta method. We find that MC results 
(both in the MCE and CE) and the MF results are in excellent agreement. (This happens although the geometry for 
the MC calculation is cubic while it is spherical for the MF). 

The clumping phase transition takes place when Dc(tj r ) vanishes at 77^ = rj R . At such point the expansion in X/N 
breaks down since the correction terms become large in eqs. (14-15]). MF applies when N\rj R — n R \ ^> 1. Since, 
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eq.(|l^) correctly suggests that PV/[NT] becomes large and negative for 77^ f 77^ in agreement with the MC results. 
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FIG. 1: f MF (v R ) =PV/[NT] as a function of r) R in the MF approximation [eq.(|l5|)], Do(ri R ) for the canonical ensemble (CE) 
and D% IC for the microcanonical ensemble (MCE). 



From eqs . (JTl|-|T5|) we obtain the following behaviour near the point C in fig.l 

f MF (n*) ^ 1/3 + 0.213738 . . . y/tjg - r, R + Ofjfe - v R ) 
{c v ) M F 0.80714... (^-^)"V2_ 0.19924... + 0{JrfT^) 
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As noticed before, the CE only describes the region between the points IG (ideal gas point, rj = 0) and C in fig.l. 
The MCE goes beyond the point C (till the point MC) with the physical magnitudes described by the second sheet 
of the square root in eqs.(|l7|). We have between C and MC 

f M F(v R ) ^ 1/3 - 0.213738 . . . sffe - if + 0( V % - if) 

(c v ) MF -0.80714 ...(rig- ii R y 1/2 - 0.19924 . . . + 0{yjr,g - r, R ) 

cy is negative in all the interval from C to MC and vanishes at the point MC. 

A functional integral representation for w(£,N) in the MCE follows by inverse Laplace transform of e*™^ ^ [0. 
The saddle point in the MCE is Ns(r) R ) as in in the CE but the corrections are different. The S-wave determinant 
results, D S MC = 6f MF {r) R ) 2 - (11/2 - V R )Im F {v r ) + 1/2. It is positive from IG to MC where it vanishes (fig.l). 

Speed of sound and compressibility. 

The isothermal compressibility and the speed of sound for long wavelengths [|| follow from the equation of state (||) 
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K F is positive for < r) R < rj R = 2.4345... where K F diverges, it is negative for r/ R < r\ R < rj R . The singularity of K F 
before but near the point C appears as a preliminary signal of the phase transition and perhaps 77,^ is the transition 
point T seen with the Monte Carlo simulations (sec fig. 2). K F becomes positive between C and MC. 

The speed of sound squared v^/T(rj R ), is positive and decreasing in the whole interval between / and C. At the 
point C it takes the value v 2 /T(r/ R ) = 11/18. Then, v^/T(r] R ) decreases between C and MC becoming negative at 
r) R = 2.14674... v 2 < indicates an instability announcing the MC critical point at i] R c — 2.03085... 
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FIG. 2: f{rj) = PV/[NT] as a function of -q by Monte Carlo for the MCE and CE (N = 2000). 



Particle Distribution. 

In the dilute regime r] R <C 1 the gas density is uniform, as expected. We find that these mass distributions approxi- 
mately follow the power law 



M(r)~Cr 
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where D slowly decreases with X(r] R ) as depicted in Table 2 from the value D = 3 for the ideal gas (77 = 0) till 
D = 0.98 in the extreme limit of the MC point. 

The gravitational potential at the point r in the MF follows from eq. (|Io| ) to be 

U{r) = --[4>{x)-a] . 
m 

The local pressure in the gas takes then the form p(r) = T p{r). We have thus derived the equation of state for the 
self-gravitating gas. It turns to be locally the ideal gas equation. But, the self-gravitating gas is inhomogeneous, 
the pressure at the surface of a given volume is not equal to the temperature times the average density of particles in 
the volume. In particular, for the whole volume: PV/[NT] — f(rj R ) ^ 1 (except for 77 = 0). 
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TABLE 1. The Fractal Dimension D and the proportionality coefficient C as a function of rj R from a fit to the 
mean field results according to M(r) ~ C r D . 

We have thus shown that the thermodynamics of the self-gravitating gas in the N — > 00, L — » 00 limit with N/L 
fixed can be expressed in terms of a single function f(rj). Besides computing/ (77) numerically we obtained its Riemann 
sheet structure and its behaviour near the critical point analytically [eq. (|17j)] . 
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